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ABSTRACT 

Direct  integration  leads  to  the  potential  gradient  of  a  rectangular  plate.  Lines  with 
equal  spacing  form  a  grid  for  integration  over  the  infinite  plane  A  source  distribution 
is  approximated  bv  pattern  functions  which  are  centered  over  each  grid  point  The 
pattern  functions  are  based  upon  the  sine  quotient  function  Fourier  analysis  replaces 
integration  in  physical  space  with  integration  in  wave  number  space  The  range  of 
integration  is  limited  to  lines  of  finite  length  in  wave  number  space  The  computation 
of  potential  gradient  is  provided  by  subroutines 


INTRODUCTION 


The  flow  around  a  ship  is  partly  a  radial  flow  from  a  source  distribution  and  partly 
a  circular  flow  around  a  vortex  distribution  In  either  case  the  velocity  is  the  integral 
of  the  product  of  an  inverse  square  and  the  density  of  sources  or  vortices  The  inverse 
square  contributes  a  spike  to  the  integrands  and  numerical  integration  is  not  feasible 
A  technique  for  integration  through  the  spike  may  be  found  from  an  analysis  of  the 
analogous  problem  in  the  flow  over  a  plane 

The  potential  and  the  potential  gradient  of  a  source  or  charge  which  is  distributed 
over  a  plane  have  applications  in  hydrodynamics,  electromagnetism,  and  gravitation 
The  field  of  an  irregular  three-dimensional  source  distribution  can  be  constructed 
with  source  distributions  over  each  of  a  stack  of  planes.  The  potential  and  the  potential 
gradient  at  a  field  point  near  a  plane  can  be  computed  by  two  methods. 

In  a  direct  method,  the  potential  at  a  field  point  is  the  integral  over  the  plane  of 
the  product  of  inverse  distance  and  the  source  density,  while  the  potential  gradient 
is  the  integral  over  the  plane  of  the  product  of  an  inverse  square  of  distance  and  the 
source  density  Wherever  there  is  a  discontinuity  in  source  density  the  potential 
gradient  is  infinite 

In  a  Fourier  method,  the  potential  is  expressed  as  a  solution  of  Laplace's  equation, 
and  the  solution  is  adjusted  to  fit  the  boundary  conditions  The  solution  is  the  sum 
of  products  of  complex  exponential  functions,  and  the  source  distribution  is 
approximated  as  a  Fourier  series  A  Fourier  transform  replaces  integration  in  physical 
space  with  integration  in  wave  number  space  The  time  of  computation  depends  upon 
the  distance  to  the  field  point  There  are  two  methods  for  the  integration  in  wave 
number  space  In  the  first  method,  Gauss  integration  multipliers  are  applied  in  a 
roving  interval  The  time  increases  with  distance  In  the  second  method,  integration 
by  parts  leads  to  recurrence  equations  The  time  is  independent  of  distance 


DIRECT  INTEGRATION 


Let  x.y.z  be  Cartesian  coordinates  in  a  right-handed  coordinate  system  with  x  and  y 
in  the  plane  and  with  z  above  the  plane  Let  a(x.y)  be  the  source  density  at  a  point 
with  coordinates  x,  y  Then  the  potential  is  given  by  the  equation 


y.  z)  = 


o(u,  v) 

\f(x  -  u)s  +  (y  -  v)z  +  zz 


dudv 


(1) 


where  x.y.z  are  coordinates  of  a  field  point  and  u.v  are  coordinates  of  a  source  point 


RECTANGULAR  PLATE 


Let  a  rectangular  plate  be  placed  on  the  plane  Let  the  length  of  the  plate  be  2a. 
and  let  the  breadth  of  the  plate  be  26  Let  the  origin  be  at  the  center  of  the  plate 
with  the  x  axis  in  the  direction  of  length  and  with  the  y-axis  in  the  direction  of 
breadth 

The  potential  of  a  plate  with  unit  source  density  is  given  by  the  equation 


<fi(x.  y.  z ) 


'*b  y  /■*•-*  _  dudv 

1  -6  y  1  ■  o  -  I  V  11  4  r 


(2) 


Initial  integration  leads  to  inverse  hyperbolic  functions,  then  final  integration  is 
completed  with  integrations  by  parts 


V  V 


The  potential  is  given  by  the  equation 

/  x  u-i  (b  ~y)  ,  /  \  .  i_  - 1  (b  +  y) 

<f  =  (a  -  x)sinh  '—============  + (a  -  z)sinh  -  ==  == 

v  (a  -  x)2  +  zz  v(a  x)2  4  2  2 


(6  -  y) 

4  (a  4  x)sinh  1  — ~~ — —  ... .  4  (a  4  x)sinh' 

v  (a  4  x)2  4  zz 


i  (b  +  y) _ 

V(a  4  x)2  4  22 


,  /t  x  u-,  (a  -  x)  ,  x  ■  l  —  i  (a  +  x) 

+  (6-!/,s,nh  7 E  -  V)»  V?  + (b ' a>smh 

+  (6  +  j/)sinh*‘  -■ - 2 - +  (6  +  y)sinh~'  y-  i  ^ 

V(bT^)z  +  z2  y  v^(6  4y)24  2a 

- 1  (a  -  x)(6  -  y)  (a  -  x)(6  4  y) 

-  z  tan  —  —  - -  -  z  tan  — j— - ^rr . 

:v(a  -  x)2  4  (6  -  y)2  4  z2  zV(a  -  x)2  4  (6  4  y)2  4  z2 


,  (a  4  x)(6  -  y)  _ 

-  x  tan  — 7===^  -  x  tan  — 

2\/(a  4  x)2  4(6-  y)2  4  22  2 


_ (a  4  x)(6  4  y) _ 

^(a  -I  x)2  4(6+  y)2  4  22 


Partial  difTerentiation  and  cancellation  of  terms  gives  the  components  of  the  gradient 
of  the  potential 

Differentiation  with  respect  to  x  leads  to  the  equation 


av>  u  ,  («>  -  y)  ■  ,  (6  +  y) 

-  -  =  sinh  — r- - ===— -  +  sinh  — - 

ajr  v^cr  -  x)2  +  z2  \/{a  -  x)2  4  22 


,  ,  (by)  ,  (6  4  y) 

sinh  '——--=====  -sinh  1  —7=^—  - - 


\/(a  +  x)2  +  2 2 


V(a  4  x)2  4  2 2 


differentiation  with  respect  to  y  leads  to  the  equation 

df  ,  (a  -  x)  (a  4  x) 

-  --  “  sinh  —  ==■ ■__=  ■ -  =  +  sinh  —  ,  - 

ay  v(6  y)2  4  22  V(6-v)2  + 


y)2  +  2 2 


u  ,  (a  -  x)  (a  4  x) 

-Sinh  -  .  ~  SUlH  '—7= . .  . - 

v(6-ty)2  +  22  v(6  4  y)2  +  2 2 

and  differentiation  with  respect  to  z  leads  to  the  equation 


(a  -  x)(6  -  y) 

-  tan  1  — —  — - - - —  - - r- 

z-J {a  x)2  t  (6  -  y)2  +  22 

(a  4  x)(6  -  y) 

4  tan  1 — v  --  — - .-=== 

2v(a  4  x)2  +  (6  -  y)2  4  22 


4  tan 


, _ (a  -  x)(b  4  y) 

;V(a  -  x)2  4  (6  4  y)2  4  2 2 


(a  4  x)(6  4-  y) 
4  tan  1 — — - .  - 


zv(a  4  x)2  4  (6  +  y)2  4-  2* 


These  equations  were  published  in  a  previous  report1  They  are  used  m  the  following 
subroutines 
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SUBROUTINE  RPLTP  (AA,  AB.  AX,  AY,  AZ.  FP) 

FORTRAN  SUBROUTINE  FOR  POTENTIAL  OF  RECTANGULAR  PLATE 


The  half-length  a  of  the  plate  is  given  in  argument  AA,  the  half-breadth  6  of  the 
plate  is  given  in  argument  AB,  and  the  Cartesian  coordinates  x,  y,  z  of  the  field  point 
are  given  in  arguments  AX,  AY,  AZ  The  potential  of  the  plate  is  stored  in  function  FP. 


SUBROUTINE  RPLTG  (AA,  AB,  AX,  AY,  AZ,  FX,  FY,  FZ) 

FORTRAN  SUBROUTINE  FOR  POTENTIAL  GRADIENT  OF  RECTANGULAR  PLATE 
********************************************************************************************* 


The  half-length  a  of  the  plate  is  given  in  argument  AA,  the  half-breadth  6  of  the 
plate  is  given  in  argument  AB.  and  the  Cartesian  coordinates  x,  y,  z  of  the  field  point 
are  given  in  arguments  AX,  AY,  AZ  The  components  of  the  gradient  are  stored  in  functions 

FX, FY, FZ 


NONUNIFORM  PLATE 


A  major  achievement  was  the  integration  when  the  source  density  was  expressed 
as  power  polynomials  in  the  surface  coordinates.  The  computation  of  potential  required 
a  triplex  of  recurrence  relations.  The  computation  is  provided  by  the  subroutine  RPLTM 
The  analysis  has  been  presented  in  a  previous  report1  That  RP1.TG  and  RPLTM  give  the 
same  gradients  has  been  verified  by  computation 


FOURIER  ANALYSIS 

The  Fourier  transform  for  a  function  on  a  plane  is  given  by  the  equations 


A(a,  p)  =  J  J  F(x,  y)e~'lax*tv)  dx  dy 


(7) 


F(x.  y)  =  J  J  Aia.  P)  e'{az*fv)  da  dp 


(8) 


where  x.y  are  Cartesian  coordinates  in  physical  space,  a.p  are  Cartesian  coordinates 
in  wave  number  space.  F(x,  y)  is  a  function  in  physical  space,  and  A(a.  p)  is  the  Fourier 
amplitude  in  wave  number  space  A  potential  ip{x,  y.  z)  of  a  source  distribution  on  the 
plane  can  be  constructed  with  the  equation 


<fi(x,  y.  z)  =  J  J  A(a.p)e  ■  da 


dp 


(9) 


which  gives  a  solution  of  Laplace's  equation  wherever  2*0  The  vertical  component 
of  the  gradient  at  z  =  0  is  given  by  the  equation 


difi 


dz 


\  ~p~z  A(a.p)e'{al'fv>  da  d<s 


(z  -  0)  (10) 


where  the  sign  is  +  if  2  >  0  and  the  sign  is  -  if  z  <  0  In  accordance  with  the  Gauss 
theorem  the  difference  in  d<p/dz  on  opposite  sides  of  the  plane  is  4no(x,  y)  where 
a(x,  y)  is  the  surface  source  density  at  the  plane  Thus  the  Fourier  amplitude  for  an 
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arbitrary  distribution  of  source  density  is  given  by  the  equation 

1 


A(a.p)  =  J  f  "(X'y)e-ilai+'v)dxdy 


(11) 


2nV  a2  +  pz 

The  range  of  integration  of  density  is  the  infinite  plane 

RECTANGULAR  PLATE 

Let  a  rectangular  plate  be  placed  on  the  plane.  Let  the  length  of  the  plate  be  2a, 
and  let  the  breadth  of  the  plate  be  26  Let  the  origia  be  at  the  center  of  the  plate 
with  the  t  axis  in  the  direction  of  length  and  with  the  y-axis  in  the  direction  of 
breadth  Let  u.  v  be  the  coordinates  of  a  source  point  on  the  plate  and  let  x.y.z  be 
the  coordinates  of  a  field  point  above  the  plate. 

The  potential  of  a  plate  with  unit  source  density  is  given  by  the  equation 

-  •<ja*<-02lll  +  i|a(x-u>+0(y-v)| 


y.  Z)  = 


1 


JJJ7 


2  n  J  J  J  J  V^Tp* 

Completion  of  the  integration  leads  to  the  equation 


da  dp  du  dv 


2  f  f  sin(aa)  s 

y.  z)  -  -  — - - 

n  i  J  a 


sin(6/3)  e-x/a2^2|ll+t<ax4/,v) 


P 


\/a2  +  pz 


da  dp 


(12) 


(13) 


What  is  interesting  about  this  equation  is  that  the  distribution  of  source  density  over 
a  finite  plate  introduces  the  factors  sin(aa)/a  and  sin(bp)/p  into  the  Fourier  amplitude 
The  factors  concentrate  the  spectrum  of  the  Fourier  integral  into  bands  which  are 
centered  with  respect  to  the  line  spectrum  of  a  trigonometric  polynomial 

INTERPOLATION 

The  coefficient  function  for  Lagrange  interpolation  in  an  infinite  set  of  equally 
spaced  data  is  just  the  infinite  product  for  the  sine  quotient  function  in  the  equation 


^  =  n  (,  -  ^ 

XI  \  kz n! 


(14) 


It  is  equal  to  unity  at  the  origin  and  is  equal  to  zero  wherever  else  u  is  a  multiple 
of  tt.  It  is  analytic  and  diminishes  with  distance  from  the  origin 

An  arbitrary  function  /(it)  is  expressed  in  terms  of  coefficient  functions  with  centers 
at  kn  by  the  equation 


,,  ,  ^  ,  sin(u  -  kv) 

f{u)  =  £  f(ut) - - - 

u  -  k TT 


(15) 


where  f(uk)  is  the  discrete  value  of  f(u)  at  the  kth  node 

The  product  of  functions  which  are  centered  at  kn  and  mn  is  given  by  the  function 


sin(u  k-n)  sin(u  -  mn) 


(16) 


(u  ~  fcn)  (u  -  mn) 

Resolution  into  partial  fractions  and  application  of  the  addition  theorem  for 
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trigonometric  functions  leads  to  the  equation 
f+“  sin(u-Jt7T)  sin(u-mTT)  ^  _  ( —  1  )*+m T  T+0°  sina(u-Jfc7r)  r4 

J-oo  (u-ki t)  (u-tnn)  (A:-m)n J_«„  u-kn 


sina(u-mn) 


du  (17) 


Inasmuch  as  the  integrands  are  odd  their  integrals  are  zero  and  the  coefficient 
functions  are  orthogonal  with  respect  to  integration  They  may  be  used  as  the  basis 
for  approximations  by  least  squares 

An  important  property  is  expressed  by  the  equation 


-n  if  a  <  0 


0  if  a  =  0 


+  n  if  a  >  0 


where  a  is  a  constant  of  arbitrary  magnitude  Integration  by  parts  leads  to  the  equation 


f+”  sinau  f  *"  sin  2u 

- r  du  =  - du  = 

J—  ua  u 


which  gives  the  normalization  factor  for  least  squares  approximation 


SOURCE  PATTERN 


A  universal  source  pattern  for  the  infinite  plane  is  given  by  the  equation 


F(x,  y)  = 


'MtM 

mm 


The  function  F(x,y)  is  unity  at  the  origin  and  is  zero  on  all  equally  spaced  nodal  lines 
which  do  not  pass  through  the  origin  The  spacing  between  nodal  lines  is  2a  in  the 
i-direction  and  is  26  in  the  y~direction  The  function  is  analytic  and  diminishes  with 
distance  from  the  origin  Application  of  the  Fourier  transform  and  use  of  the  Euler 
theorem  introduces  the  difference  of  two  cosines  and  introduces  the  sum  of  two  sines 
into  the  integrand  of  the  Fourier  amplitude.  The  difference  between  the  cosines  is 
zero  at  the  origin  and  is  an  even  function  otherwise.  The  difference  between  the  cosines 
is  eliminated  by  integration  There  remains  the  sum  of  sines  as  expressed  by  the 
equation 


=  J_  f  f  sin(^-a)z  +  sin(^,  +  q)x  sin(^-fi)y  +  sin(^  +  ff)y 

*-wJJ  ^  m 


m 


The  sines  cancel  each  other  during  integration  wherever  a.fi  are  outside  a  rectangle 
in  wave  number  space  The  interior  of  the  rectangle  is  located  where  a,  (}  satisfy  the 
inequalities 


7T  TT 

-  iqs  +  - 

2a  2a 


5 


Vi 


Thus  the  function  F(x,  y)  is  given  by  the  equation 


F(x.y) 


_  ab  rs  r 

J  '  J. 


dcx  dp  =  o(x.  y) 


and  the  potential  x .  y.  z)  is  given  by  the  equation 


rt-r.  y.  =) 


2  a  b 
rr 


Zb 


z  *\lax*0y) 


da  dfi 


(23) 


(24) 


Further  integration  is  possible  after  Cartesian  coordinates  n  fi  are  replaced  bv  polar 
coordinates  k.  6  The  potential  is  expressed  by  the  equation 


y  U-  y  z) 


2ab 

7T 


*  z  *  ix  *  rcosfl  •  y  #  i  n  0 1 


(25) 


Then  integration  with  respect  to  k  replaces  surface  integration  throughout  the  interior 
of  the  rectangle  with  line  integration  along  the  perimeter  of  ttie  rectangle 

The  three  components  of  the  gradient  are  given  by  the  same  formulations  except 
for  a  weighting  function  u  Thus  for  the  components 

<3y  _  d*r 

d-r  dy 

the  values  of  u  ai  > 

i  cos  f‘  -  i  sm  6 

Then  each  component  is  the  sum 

.4  •  B 

of  two  integrals  where  4  is  the  integral  along  vertical  sides  of  the  rectangle  and  B  is 
the  integral  along  horizontal  sides  of  the  rectangle 

For  the  vertical  sides  the  variables  are  related  in  accordance  with  the  equations 

I  tan  6  (39) 


dy 

dz 


(26) 


(27) 


(28) 


sm  0 


COS  6  -  ~  -----  =:  — . 

\  1  -  t~ 


dS 


dt 


(30) 


o  «.  cos  6  -  (31) 

2  a 

Let  a  parameter  3  tie  defined  by  the  equation 


(5  !  :  \  !  •  f;  i ( x  -  y f I i  (32) 

da 

Then  the  integral  4  is  given  bv  the  equation 

b  i"s  I  (  1  •  3)e  6 

4  ~  I  .  n  dt  (33) 

a  a  o“ 


for  integration  along  'mt h  vertical  sides 


6 


For  the  horizontal  sides  the  variables  are  related  in  accordance  with  the  equations 

t  -  cot  0  (34) 


sin  6  = 


1 


Vl  H  tZ 


cos  6  = 
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v/TTT1 


^  =  ac  sin  0  =  - 

Let  the  parameter  <5  be  defined  by  the  equation 


dec¬ 


eit 


1  +  f* 


(35) 

(36) 


6  =  —  ||z|n/T+  tz  -  i(xt  +  y)\ 

co 


Then  the  integral  B  is  given  by  the  equation 


B 


a  C*a  1  -  (1  +  <5)e~a 


iu  dt 


(37) 


(38) 


for  integration  along  both  horizontal  sides. 

Wherever  j<5 1  S  1  the  integrand  is  replaced  by  the  infinite  series  in  the  equation 


-£(*+■>  (~6)t 


k  =  0 


(*  +  2)! 


(39) 


Summation  of  this  series  is  continued  until  there  is  no  change  in  the  sum 


GAUSS  INTEGRATION 


The  integration  is  completed  with  the  aid  of  16-point  Gaussian  integration 
multipliers  Each  interval  of  integration  is  subdivided  into  subintervals  and  the  Gaussian 
integration  is  applied  progressively  to  successive  subintervals  The  numbers  of 
subintervals  for  the  integration  of  A  and  B  are  the  integral  parts  of  the  expressions 


2  + 


lyl 
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(40) 


Then  the  Gaussian  integration  is  accurate  for  any  x  or  y  For  large  values  of  x,  y  the 
time  for  computation  is  nearly  proportional  to  the  sum  of  \x\  and  |y| 

Formulae  for  the  Taylor  series  expansion  of  the  source  pattern  at  nodal  points  in 
a  finite  grid  were  used  with  subroutine  RPLTM  to  obtain  gradients  at  several  grid  sizes 
Extrapolation  to  infinite  grid  size  confirmed  the  line  integration  in  wave  number  space 
The  components  of  the  gradient  of  potential  are  given  by  the  following  subroutine 

SUBROUTINE  rptng  (aa,  ab,  ax,  ay,  a z,  px.  fy,  fz) 

FORTRAN  SUBROUTINE  FOR  POTENTIAL  GRADIENT  OF  RECTANGULAR  PATTERN 


The  half-interval  a  of  the  pattern  is  given  in  argument  AA,  the  half-interval  6  of  the 
pattern  is  given  in  argument  AB,  and  the  Cartesian  coordinates  x,  y,  z  of  the  field  point 
are  given  in  arguments  AX,  AY,  AZ  The  components  of  the  gradient  are  computed  by 
Gauss  integration  and  are  stored  in  the  functions  FX,  FY,  FZ 
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INTEGRATION  BY  PARTS 


For  the  integration  by  parts,  the  quadratic  formula  is  used  to  solve  the  equations 
for  6  to  give  t  as  a  function  of  6 

For  the  vertical  sides  of  the  rectangle  the  variable  t  is  given  by  the  equation 


iy(^  +  **)  *  +  ix)2  -  (y2  +  z2) 

V*  +  za 


(41) 


Substitution  and  cancellation  lead  to  the  equation 


\fi  +  t2 


IfK^jr  +  ix)  ±  -*•  ix)2  -  ( y z  +  zz) 


yr  +  z 


(42) 


and  differentiation  leads  to  the  equation 

dt  _  H  \iyy^*ix)t-(yz^zr)  ±\z\(2-f  +  ix) 
d6  (y2  +  z2)  [  +  ix)2  -  (y2  +  zz) 


(43) 


The  ±  sign  in  the  equations  is  determined  by  substitution  of  the  limits  of  integration 
for  t  The  sign  is  +  for  positive  t  and  is  -  for  negative  t 

Substitution  in  the  equation  for  A  replaces  integration  with  respect  to  t  with 
integration  with  respect  to  <5  The  natural  path  of  integration  in  the  f-plane  is  the 
real  axis,  whereas  the  natural  path  of  integration  in  the  6-plane  is  an  hyperbola  with 
vertex  where  6  is  given  by  the  equation 


6  = 


TT 

2a 


itzl  -  tx| 


(44) 


and  with  axis  in  the  direction  of  the  positive  real  axis  There  is  a  singularity  on  each 
side  of  the  imaginary  axis  Near  the  singularity  on  the  positive  side  6  is  given  by  the 
equation 


<5  =  5-  1  v  y 

2a 


/V  + 


ix  +  e 


(45) 


Then  the  limiting  values  for  functions  at  the  singularity  are  given  by  the  equations 


v'y2  *  z2 


/  2  2 
v  y  +  zc 


and  by  the  equation 


\  < 


dt 

d6 


z  v  , 


( y 2  *  z2)* 


(46) 


(t  -  0)  (47) 


The  hyperbolic  path  in  the  6-plane  is  deformed  into  two  straight  lines  which  run  from 
the  positive  singularity  to  the  limits  of  integration  where  /tab 

For  the  horizontal  sides  of  the  rectangle  the  variable  t  is  given  bv  the  equation 
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i*( 


Zbd 


xy)  t  s  n  (2?  *  <?/>  = 
x 2  *  z2 


z2) 


(48) 
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Substitution  and  cancellation  lead  to  the  equation 

IzKNr  +  w)  ±  ix\ /(H*  +  w)z  -  (*z  +  22) 


and  differentiation  leads  to  the  equation 

“  +  iy)2  -  (x2  +  22)  ±|*|( 


d6  (x2  +  22)  L  V(^2  +  ly )2  -  (x2  +  2 2) 

The  ±  sign  in  the  equations  is  determined  by  substitution  of  the  limits  of  integration 
for  t  The  sign  is  +  for  positive  t  and  is  -  for  negative  t 

Substitution  in  the  equation  for  B  replaces  integration  with  respect  to  t  with 
integration  with  respect  to  6  The  natural  path  of  integration  in  the  f-plane  is  the 
real  axis,  whereas  the  natural  path  of  integration  in  the  d-plane  is  an  hyperbola  with 
vertex  where  6  is  given  by  the  equation 

<5  =  fill  -  »V|  (51) 

cO 

and  with  axis  in  the  direction  of  the  positive  real  axis  There  is  a  singularity  on  each 
side  of  the  imaginary  axis  Near  the  singularity  on  the  positive  side  6  is  given  by  the 
equation 

<5  =  ~  +  22  -  iyi  +  c  (52) 

cb 

Then  the  limiting  values  for  functions  at  the  singularity  are  given  by  the  equations 


'X2  +  2 2 


'1  +  f2  = 


t X 2  +  22 


and  by  the  equation 


/—  at 
vc  —  -* 


(x2+  22)4 


(e  -  0)  (54) 


The  hyperbolic  path  in  the  (5-plane  is  deformed  into  two  straight  lines  which  run  from 
the  positive  singularity  to  the  limits  of  integration  where  t  =  ±  6/a 
The  weighted  derivative 


is  expressed  as  a  power  series  expansion  by  11 -point  Lagrange  interpolation  Let  6 
be  given  by  the  equation 

6  =  t)  +  c  (56) 

where  rj  is  a  center  of  expansion  and  c  i  the  argument  of  expansion  Inasmuch  as 
the  paths  of  integration  are  straight  lines  in  the  f-plane,  all  values  of  £  in  the  range 
of  interpolation  are  given  by  the  equation 

£  =  U(M  (-1  SuS  +1)  (57) 

where  is  the  upper  limit  of  interpolation  and  u  is  the  variable  of  interpolation  The 
limit  t*  can  be  cancelled  out  of  the  rational  expression  for  the  complex  Lagrange 


fit 


coefficient  and  the  derivative  can  be  expanded  as  a  power  series  in  u.  In  a  preliminary 
computation  discrete  values  for  u  with  Chebyshev  spacing  and  discrete  coefficients 
were  computed  for  use  in  a  subroutine.  Thus  the  derivatives  are  expressed  as  power 
polynomials  in  the  ratio 

u  =  —  (58) 

e* 

while  division  by  eM  is  deferred  until  the  integration  by  parts. 

There  are  three  stages  of  expansion.  In  the  first  stage,  17  coincides  with  the  singularity 
and  derivatives  are  expanded  as  power  polynomials  in.£,/J  The  range  of  interpolation 
is  extended  into  a  negative  range  in  order  to  obtain  data  for  central  interpolation. 
The  extension  in  the  e1/z-plane  corresponds  to  a  fold  in  the  e -plane.  The  interval  of 
integration  is  the  upper  half  of  the  range  of  interpolation.  In  the  later  stages  of 
expansion.  r\  is  moved  further  out  along  the  path  of  integration  and  the  derivatives 
are  expanded  as  power  polynomials  in  e  The  interval  of  integration  is  the  full  range 
of  interpolation 

There  are  tight  limitations  on  the  range  of  interpolation  The  radius  of  expansion 
near  the  first  singularity  must  be  a  small  fraction  of  the  distance  from  the  first 
singularity  to  the  second  singularity  The  radius  of  expansion  near  the  second  center 
of  expansion  must  be  a  small  fraction  of  the  distance  from  the  center  of  expansion 
to  either  singularity 

The  integration  of  each  series  in  powers  of  £1/z  is  completed  with  a  descending 
recurrence  Required  for  the  integration  are  integrals  which  are  generated  by  the 
recurrence 
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The  iteration  is  started  with  an  initial  approximation  for  m  =  32.  Completion  of 
integration  is  given  by  integrals  which  are  obtained  by  the  recurrence 
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where  m  is  alternately  half  an  odd  integer  and  half  an  even  integer.  The  iteration  is 
started  with  an  initial  approximation  for  m  =  32  Stability  of  the  recurrence  is  achieved 
by  a  limitation  of  |e|  to  ^g|7j|. 

The  integration  of  each  series  in  powers  of  e  is  cycled  in  either  direction  Required 
for  the  integration  are  integrals  which  are  generated  by  the  ascending  recurrence 
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The  iteration  is  started  with  the  special  integral  in  the  equation 

f  £  e*a  d(  =  e~v  -  (1  +  £)e~a  (62) 

Jo 

and  the  iteration  is  continued  until  m  =  64  If  |e|  §  17  then  the  iteration  is  cycled  in 
descending  order  with  an  initial  approximation  Completion  of  integration  is  given  by 
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integrals  which  are  generated  by  the  ascending  recurrence 
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where  m  is  a  sequence  of  integers  The  iteration  is  started  with  the  aid  of  the  special 
integrals  in  the  equations 


£' '  -  < 1  *  "w“  -  -  *  ♦  *  *  <*  ♦ «»-  -  -  r  ^rmrw.  (66> 

and  the  iteration  is  continued  for  1 1  cycles  If  |c|  £  ||ry|  the  iteration  is  cycled  in 
descending  order  with  an  initial  approximation 

In  the  actual  program  e  is  replaced  by  t/c*.  which  gives  integrals  ready  for  use 
with  the  Lagrange  coefficients,  and  scales  the  integrals  so  that  their  exponents  do  not 
exceed  the  tight  limitation  on  the  exponents  in  IBM  computers.  The  recurrence  equations 
can  be  verified  easily  by  differentiation 

Although  the  integration  by  parts  is  accurate  everywhere  that  \z\  is  not  small,  the 
integration  by  parts  breaks  down  when  yz  +  z2  is  zero  along  the  vertical  sides  and 
breaks  down  when  z2  +  z2  is  zero  along  the  horizontal  sides  Under  such  circumstances 
the  program  falls  back  upon  Gauss  integration 

For  large  values  of  z,  y  the  time  for  computation  is  independent  of  |x|  or  |y|.  The 
components  of  the  gradient  of  the  potential  are  given  by  the  following  subroutine. 

SUBROUTINE  rptnx  (aa,  ab.  ax,  ay.  az,  ex,  fy,  FZ) 


y  (-6)'" 
*.0  (*  +  2 )• 


1  +  e~4  +  y  +  log  <5  -  Ei(-6)  =  £ 


k%  (k  +  2)zk\ 


FORTRAN  SUBROUTINE  FOR  POTENTIAL  GRADIENT  OF  RECTANGULAR  PATTERN 


The  half-interval  a  of  the  pattern  is  given  in  argument  AA,  the  half-interval  6  of  the 
pattern  is  given  in  argument  AB,  and  the  Cartesian  coordinates  z,  y.  z  of  the  field  point 
are  given  in  arguments  AX.AY.AZ  The  components  of  the  gradient  are  computed  with 
integration  by  parts  and  are  stored  in  functions  FX,  fy,  FZ. 


GRADIENT  OVER  GRID 

Exploratory  computations  have  given  the  gradient  of  potential  on  a  vertical  line 
through  each  grid  point  of  the  source  pattern  The  gradient  has  a  direction  which 
passes  through  the  center  and  has  a  magnitude  which  varies  with  height  Far  from 
the  plane  the  magnitude  approaches  the  full  value  for  the  inverse  square  law.  while 
near  the  plane  the  magnitude  is  diminished  by  shielding  by  the  local  structure  of  the 
source  pattern  As  distance  diminishes  toward  the  center  the  inverse  square  of  distance 
approaches  °°,  while  the  gradient  approaches  2tt 


V/. 


The  potential  is  given  by  the  equation 
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da  dp 
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Combinations  of  the  partial  derivatives  replace  the  integrand  with  an  exact  differential 
which  can  be  integrated  at  once  Results  of  integration  are  given  by  the  equation 
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The  integral  is  zero  if  x  is  a  multiple  of  2a.  Results  of  integration  are  given  by  the 
equation 
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The  integral  is  zero  if  y  is  a  multiple  of  26  Thus  the  potential  gradient  has  a  radial 
direction  over  grid  points  where  x.  y  are  multiples  of  2a,  26  The  integration  of  the 
inverse  square  over  the  surface  of  the  pattern  is  just  the  gradient  of  a  point  source 
when  the  field  point  is  over  a  grid  point 


RECTILINEAR  LINE 

When  the  plane  contracts  to  a  line,  the  potential  and  the  potential  gradient  become 
symmetric  with  respect  to  the  line  The  position  of  a  field  point  is  located  by  the 
cylindrical  polar  coordinates  r,  0,  z,  where  r  is  the  distance  from  the  line,  <f>  is  the 
azimuth  angle,  and  z  is  the  zenith  distance  along  the  line.  The  Fourier  transform  for 
a  function  on  the  line  is  given  by  the  equations 


A(ac)  =  ^  J  /(z)  e'1**  dz 

(70) 

p  +  co 

/(z)  =  A(k)  e1**  dK 

U  —  oo 

(71) 

where  z  is  the  coordinate  in  physical  space  and  ac  is  the  coordinate  in  wave  number 
space 

A  universal  source  pattern  for  the  infinite  line  is  defined  by  the  equation 


where  the  spacing  between  nodal  points  is  2c  Application  of  the  Fourier  transform 
and  use  of  the  Euler  theorem  introduces  the  difference  of  two  cosines  and  the  sum 
of  two  sines  into  the  integrand  of  the  Fourier  amplitude  The  difference  of  cosines  is 
an  even  function  of  z,  and  the  sum  of  sines  is  an  odd  function  of  z  Both  functions 
are  divided  by  the  odd  function  z  The  cosines  are  cancelled  by  the  integration  There 

12 


&&&&&& 


remains  the  sum  of  sines  as  expressed  by  the  equation 

sin(^  -  k)z  +  sin(£  +  k)z 
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dz 


Then  the  sines  cancel  each  other  whenever  x  is  outside  the  range 
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and  the  amplitude  is  given  by  the  equation 


(73) 


(74) 


4(n)  =  - 
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whenever  k  is  inside  the  range  Thus  the  function  f(z)  is  given  by  the  equation 
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(76) 


where  the  range  of  integration  in  wave  number  space  is  finite 

The  Laplacian  in  cylindrical  polar  coordinates  is  given  by  the  equation 

v.y, , d-*?  *  1  di .  i 

dr2  r  dr  r2  d$z  dz2 


(77) 


The  solutions  of  Laplace's  equation  are  given  in  many  texts  In  this  case  the  appropriate 
solution  contains  the  modified  Bessel  function  of  the  second  kind,  because  this  function 
diminishes  to  zero  with  increasing  argument2-3 
The  potential  is  given  by  the  equation 
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(78) 


where  A'o(xr)  is  the  Bessel  function  of  order  zero  The  Bessel  functions  satisfy  the 
equation 


dr 


A'0(xr)  x  A',  (xr) 


(79) 


where  A',(xr)  is  the  Bessel  function  of  order  one  The  Bessel  functions  satisfy  the  limits 
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A'0(xr)  -  -  log(xr) 


A  ,  (  h  r) 


k  r 
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at  small  radius  Thus  the  potential  satisfies  the  limit  equation 
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where  f(z)  is  the  linear  density  <j(z)  in  «k  eordance  with  the  Gauss  theorem 
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The  half  interval  r  of  the  pattern  is  given  in  argument  AC,  and  the  cylindrical 
coordinates  r,  z  of  the  field  point  are  given  in  arguments  AR,  AZ  The  components  of 
the  gradient  are  computed  by  Gauss  integration  and  are  stored  in  the  functions  FR,  FZ 


DISCUSSION 


The  ubiquitous  function 


is  sometimes  called  the  diffraction  function  because  it  appears  in  the  diffraction  of 
waves  at  an  aperture,  but  it  appears  in  other  applications  as  well  It  should  have  a 
name  of  its  own  It  is  proposed  that  it  be  called  the  sine  quotient  function 

The  source  density  is  a  scalar  while  the  vortex  density  is  a  vector  The  sine  quotient 
function  ran  be  used  for  the  interpolation  of  the  source  density  or  either  of  the 
components  of  vortex  density  Then  the  potential  gradient  of  the  pattern  function  is 
used  directly  with  source  density  in  the  computation  of  radial  flow,  or  in  a  vector 
produc  t  with  either  component  of  vortex  density  m  the  computation  of  circular  flow 
A  technique  for  the  computation  of  ideal  flow  over  arbitrary  bodies  has  been  invented 
by  Hess  and  Smith4’'  The  same  technique  has  been  adopted  for  computation  of  flow 
over  ships  by  Dawson  and  Dean87  In  this  tec  hmque  the  surface  of  the  body  is  subdivided 
by  a  grid,  then  each  panel  of  the  grid  is  projected  onto  a  plane  which  spans  the  four 


corners  of  the  panel.  It  is  assumed  that  the  quadrilateral  projection  has  a  uniform 
source  density.  Each  plane  quadrilateral  is  assembled  from  plane  triangles,  for  which 
the  potential  gradient  is  given  by  analytical  expressions 

In  the  case  of  an  infinite  plane,  uniform  quadrilaterals  would  be  reduced  to 
rectangular  plates  The  potential  gradient  of  a  square  plate  is  compared  with  the 
potential  gradient  of  a  source  pattern  in  the  following  table 


Potential  Gradients 

a  =  6  =  0.5  z  =  y  =  0 


z 

Inverse 

Square 

Plate 

Pattern 

0 

OO 

2n 

2rt 

0  0625 

256 

5.58063622 

5.41516931 

0  125 

64 

4.90438059 

4  68192879 

0  25 

16 

3.70918087 

3.53400628 

0  5 

4 

2.09439510 

2  09491828 

1.0 

1 

0.80543168 

0  86196571 

2  0 

0.25 

0.23543002 

0  24793473 

3.0 

0.11111111 

0.10812127 

0  11106356 

4  0 

0  0625 

0  06154089 

0  06249869 

5.0 

0  04 

0  03960461 

0  03999996 

6  0 

0.02777777 

0  02758643 

0  02777777 

7.0 

0.02040816 

0.02030466 

0  02040816 

8  0 

0.015625 

_ 

0  01556424 

0  01562500 

The  entries  in  the  table  were  computed  with  the  aid  of  RPUG,  RPTNG,  and  RPTNX  The 
table  shows  that  the  inverse  square  law  gives  accurate  gradients  for  z  >  4 

For  rectangular  plates  the  field  point  must  be  kept  away  from  the  edges,  where 
the  gradient  is  infinite,  whereas  for  pattern  functions  there  is  no  restriction  on  the 
location  of  the  field  point 


CONCLUSION 

Pattern  functions  which  are  based  upon  the  sine  quotient  function  are  useful  for 
the  computation  of  the  potential  and  the  potential  gradient  of  source  distributions 
in  the  plane  Directly  over  a  grid  point  the  potential  gradie  nt  has  a  direction  which 
passes  through  the  center  of  the  pattern,  but  has  a  magnitude  which  is  shielded  from 
the  full  inverse  square  law 


JxS: 


BIBLIOGRAPHY 

1  Computation  of  Special  Functions. 

A  V  Hershey,  Naval  Surface  Weapons  Center.  Dahlgren.  Virginia,  Report  NSWC/DL 
TR-3788  (November  1978) 

2  Modern  Analysis 

E  T  Whittaker,  and  G  N  Watson.  (Cambridge  University  Press.  Cambridge,  1952) 

3  Theory  of  Bessel  Functions. 

G  N  Watson.  (Cambridge  University  Press.  Cambridge.  1962) 

4  Calculation  of  Non- Lifting  Potential  Flow  about  Arbitrary  Three-dimensional 
Bodies 

J  L  Hess,  and  AMO  Smith.  Douglas  Aircraft  Co.,  Long  Beach.  California.  Report 
E.S  40622  (March  1962) 

5  Calculation  of  Potential  Flow  about  Arbitrary  Bodies. 

J  L  Hess,  and  AMO  Smith.  Progress  in  Aeronautical  Sciences.  8,  1  (1967) 

6  The  XYZ  Potential  Flow  Program 

C  W  Dawson,  and  J  S  Dean.  Naval  Ship  Research  and  Development  Center.  Bethesda, 
Maryland,  Report  3892  (June  1972) 

7  User's  Manual  for  Mu  XYZ  Free  Surface  Program 

B  H  Cheng,  and  J  S  Dean.  David  Taylor  Naval  Ship  Research  and  Development 
Center.  Bethesda,  Maryland,  Report  DTNSRDC -  86/029  (June  1986) 


DISTRIBUTION 

Defense  Documentation  Center 
Cameron  Station 
Alexandria,  Virginia  22314 

Library.  Code  0142 
Naval  Postgraduate  School 
Monterey.  California  93943 

Computer  Center,  Code  0141 
Naval  Postgraduate  School 
Monterey.  California  93943 

Research  Administration.  Code  012 
Naval  Postgraduate  School 
Monterey.  California  93943 

Dr  A  V  Hershev.  Code  0141 
Naval  Postgraduate  School 
Monterev.  California  93943 


